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Abstract 

The paper contains a simple approach to reconstruction in Ther- 
moacoustic and Photoacoustic Tomography. The technique works for 
any geometry of point detectors placement and for variable sound 
speed satisfying a non-trapping condition. A uniqueness of recon- 
struction result is also obtained. 
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1 Introduction 

In the past decade, one witnessed a surge in newly developing medical imag- 
ing modalities. Their designers pursue the lofty goals of increasing the image 
resolution, contrast, and safety, while reducing the costs. A new approach 
evident in some of these developments is based upon combining different 
physical types of signals in one procedure, with the hope of reducing the 
deficiencies of each individual one, and at the same time taking advantage 
of the strengths of each. The most well developed example is the Ther- 
moacoustic Tomography (TAT) (also abbreviated as TCT) [19]. Let us 
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provide a brief description of the TAT procedure (also see [HI [TSl [121 [221 [3S] 
and [in]-[lS]). 

A very short radiofrequency (RF) pulse is sent through a biological object. 
The whole object is assumed to be uniformly irradiated. Some RF energy 



is absorbed at each location inside the object. It is known (e.g., [191 141] ) 
that cancerous cells absorb several times more energy in the RF range than 
the healthy ones. Thus, if we knew the distribution of the absorption of RF 
energy function f{x), this high contrast could provide an easy identification 
of cancerous locations. However, RF waves with sufficient penetration depth 
are too long to lead to a good resolution, if used for imaging. Thus, a differ- 
ent mechanism is used for imaging f{x). Namely, energy absorption causes 
thermoelastic expansion of the tissue, which in turn leads to propagation of a 
pressure wave p{x,t). This signal is measured by transducers distributed ac- 
cording to some specifically chosen geometry around the object. Ultrasound 
imaging has rather low contrast, but in the TAT procedure, this is not a 
problem, since the contrast comes from the electromagnetic absorption, and 
ultrasound is only responsible for high resolution, which it does have. Thus, 
recovery of f{x) from the measured characteristics of the pressure p{x,t), is 
the main goal of TAT. 

The smallness of the ultrasound contrast is in fact a good thing for TAT. 
For instance, it has been mostly assumed that the sound speed is constant. 
This is the main assumption under which most known mathematical devel- 
opments have taken place. In many cases, though, it is important to take 
variable sound speed into account, and we will do so in this paper. 

The photoacoustic tomography (PAT) is almost identical to TAT. The 
difference is only in how the thermoacoustic signal is generated. In PAT, a 
laser pulse is used instead of an RF one to initiate the signal j40j. The rest 
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Figure 1: The TAT procedure. 
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of the procedure stays the same. Since from the mathematical point of view, 
there is no difference between TAT and PAT, we will be mentioning TAT 
only. 

Close attention has been paid to developing mathematical methods re- 
quired for TAT (almost solely under the constant sound speed assumption). 
The mathematical problems of TAT happen to be very interesting and in- 
volving various areas of analysis, in particular PDEs, integral geometry, mi- 
crolocal analysis, and spectral theory. Many of them, however, also are 
rather complex (see, e.e. the survey [22] and references therein)!^. After a 
substantial effort, major breakthroughs have occurred in the last couple of 
years in all relevant issues: uniqueness of reconstruction, inversion formulas 
and algorithms, range conditions, incomplete data problems, and stability. 
The hardest to come about were explicit inversion formulas, as well as some 
uniqueness of reconstruction results (albeit, for practical geometries, these 
can be considered to be resolved [3],|22]). For quite a while, only series expan- 
sions had been available for the case of transducers placed along a sphere S 
surrounding the object [211 E2]- No formulas were available for non- spherical 
surfaces 5*. The first explicit inversion formulas were obtained in [TT] for the 
spherical geometry in odd dimensions and then extended to even dimensions 
in [10]. A different (non-equivalent) set of formulas was developed for any 
dimension in [23] (see also such a formula in 3D in [32]). All these formulas 
are valid for the spherical geometry only. Another drawback of the inversions 
of [ini [m [231 132] is that if the source f{x) extends beyond the observation 
surface, the reconstruction even of its part that is inside, becomes incor- 
rect (see discussion of this phenomenon in [22]). A series expansion formula 
that involves the spectrum and eigenfunctions of the Dirichlet Laplacian was 
obtained in |25J. Unlike the formulas of [TUl \TT[ [231 [32], it applies to any 
geometry of measurement and also does not assume that the signal comes 
from the interior of the measurement surface only. It was also shown in [25] 
that the eigenfunction expansions method can be implemented to provide 
fast and accurate reconstructions. 

One needs to mention that in the case of the constant sound speed, the 
reconstruction task can be described as an integral geometry problem dealing 
with a spherical mean operator [2] 22] (see also [9l[l3l[15l[l6l[ni[23l[28l[29l[33] 
and references therein for integral geometric methods). This relation with 
integral geometry all but disappears for non-constant speed. 

^We will describe the mathematical set-up of TAT in the next section. 
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No explicit inversion formulas are known for non-spherical surfaces or for 
the case of variable sound speed. There are also some mysteries surrounding 
existing formulas. E.g., why do explicit formulas reconstruct incorrectly the 
function inside the observation surface, if the function has a part outside the 
surface? Why do inversion formulas exist, even so the general machinery of 
the so called K-operator [121 [131 [Sj suggests that one should not expect such 
formulas to appear (since one deals with a so called "inadmissible complex")? 
Some of these issues are addressed in [22]. 

The aim of this article is to develop a very simple general operator theory 
inversion formula, which applies to arbitrary geometry of the observation 
surface (i.e., the surface where transducers are placed), variable sound speed 
under a non- trapping condition, and functions not necessarily supported 
inside the observation surface. Since the formula involves functions of the 
Dirichlet Laplacian inside S, it is not easy to apply. We show how it works in 
some special cases. In particular, it leads to eigenfunction expansion methods 
that apply in a wider generality (e.g., in non- homogeneous media) than those 
of [2S]. 

One can mention that mathematical problems similar to the ones of TAT 
arise in sonar and radar research (e.g., [271 [30] ). 

The structure of the article is a s follows. The next Section [2] sets up the 
mathematical formulation of the problem. Then, in Section [3l we obtain the 
general inversion formulas. For the clarity sake, they are first derived under 
the assumption of a constant sound speed and for a domain of an arbitrary 
shape in an odd dimensional space M", so the Huygens' principle holds. Later 
on, use of the Huygens' principle is replaced by local energy decay estimates 
[381 [39]. It is also assumed that functions involved are smooth, a restriction 
that can easily be lifted after the final formulas are derived. Then an abstract 
result concerning hyperbolic equations in Hilbert spaces is derived that gen- 
eralizes the one presented in this simple example. This general formula then 
is applied to the case of reconstruction in TAT without assumption of con- 
stant sound speed. What is relevant, is not the uniformity of the background 
medium, but rather a non-trapping condition. Such a condition is absolutely 
natural, since without it one cannot expect any inverse problems of the kind 
we study to be reasonably solvable. We also establish a uniqueness of re- 
construction result. It shows that the data collected is sufficient for recovery 
of a (compactly supported) function, even if its support reaches outside the 
observation surface (albeit the inversion procedure recoveres only its interior 
part). The next Section HI contains additional remarks and discussions. This 
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is followed by an Acknowledgements section. 



2 Mathematical model of TAT 

Let the function /(x) to be reconstructed be compactly supported in M". 
Let also S" be a closed surface, at each point y of which one places a point 
detector that measures the value of the pressure p{y, t) at any moment t > 0. 
It is usually assumed (and this is crucial for the validity of the formulas of 
[ini [HI IMl 112] ) that the object (and thus the support of f{x)) is surrounded 
by S. We will not need this assumption here. 

We assume that the ultrasound speed c{x) is known (e.g., being deter- 
mined by transmission ultrasound measurements [H]). Then, the pressure 
wave p{x, t) satisfies the following problem for the standard wave equation 
[3 EH SI]: 

Ptt = c\x)A,p,t> 0,xeR^ 

pix,0) = fix), (1) 
Pt{x,0) = 

The task is to recover the initial value f{x) at t = of the solution p{x, t) 
from the measured data, which we will call g{y,t). Incorporating the data 
obtained from the measurement, the set of equations ([1]) extends to become 

'ptt = c^{x)A,p,t> 0,x G 
^ p{x,0) = f{x), 
' Ptix,0) = 

^p{y,t)=g{y,t),yeSxR+ 




Figure 2: An illustration to ([2]). 
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The problem now is to find the initial value f{x) in ([2]) from the knowledge 
of the lateral data g{y, t) (see Figure [2]). Although at a first glance one might 
think that we have insufficient data, this is in fact not true. Indeed, there is 
an additional restriction that the solution holds in the whole space, not just 
inside the cylinder S x R+. In most cases, with this additional information, 
the data is sufficient for recovery of f{x) (see [22] and references therein 
for the case of constant speed, and further in this paper for the case of a 
variable speed under a non-trapping condition). An additional sometimes 
useful comment (albeit we are not going to use it) is that p can be extended 
as an even function of time, after which it will satisfy the wave equation for 
all values of t. The same applies to the data g that can also be extended to 
an even function in t. 

We would like to notice that in the case of constant speed, due to the well 
known Poisson-Kirchhoff formulas [6l Ch. VI, Section 13.2, Formula (15)] for 
the solution of ([1]), there is close relation (essentially, an equivalence) between 
this problem and inversion of the spherical mean operator with centers of 
spheres of integration located on S (see, e.g., |31 [22] and references in these 
papers). This is why many previous considerations dealt with the spherical 
mean operator rather than the wave equation, which allowed for the well 
developed techniques (or at least hints) of integral geometry to be applied. 
We cannot relay upon such techniques, since in our case the wave speed is 
not assumed to be constant. 

The standard list of questions one would like to ask in any tomographic 
problem, includes uniqueness of reconstruction, inversion formulas and al- 
gorithms, their stability, incomplete data problems, etc. In this text, we 
will address only uniqueness of reconstruction and inversion. One can find a 
survey of known results on other issues and a bibliography in |;22|. 

3 Discussion of the problem and the main re- 
sults 

In order not to obstacle consideration, we will assume for simplicity that the 
sound speed c{x) and the initial function f{x) to be recovered are infinite 
differentiable, an assumption that can be easily significantly weakened. The 
closed observation surface S of transducers' locations is assumed to be suffi- 
ciently "nice". E.g., assuming that it consists of finitely many transversally 
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intersecting smooth pieces (e.g., a cube), would suffice both for practical ap- 
plications and for our analytic needs, albeit this condition can be weakened 
much further. Thus, we will not dwell on this issue. Notice, that unlike the 
case of all known backprojection type formulas in TAT, we do not assume 
that the function to be reconstructed is supported inside the observation 
surface S, on which the transducers are located. 

The truly significant assumptions are the following: 

1. Support of f{x) is compact. 

2. The sound speed is strictly positive c(x) > c > and such that c{x) — 1 
has compact support, i.e. c(x) = 1 for large x. 

3. Consider the Hamiltonian system in M^'l with the Hamiltonian H = 



The solutions of this system are called bicharacteristics and their pro- 
jections into are rays. 

We will assume that the non-trapping condition holds, i.e. that all 
rays (with 7^ 0) tend to infinity when t oo. 

The first assumption is very important, albeit it could be replaced by 
a sufficiently fast decay. Without a sufficiently fast decay condition, there 
is no uniqueness of reconstruction [H [3l |22] . The second condition assumes 
that the medium outside a bounded domain is homogeneous. The third 
one prohibits trapping. If trapping does occur, then it is naturally expected 
that solvability of standard inverse problems and problems like controllability 
should most probably fail. 

Let us discuss first some peculiarities of the inverse problem of recovering 
the initial value f{x) from the lateral boundary data g{y,t) in ([2]). As we 
have mentioned already, if there were nothing special about the set-up of this 
problem (i.e., if we were dealing with boundary value problem for the wave 
equation in the cylinder), such an inverse problem would not be uniquely 
solvable. Indeed, the lateral data alone is clearly insufficient. The special 
situation is that we are dealing in fact with an observability piohlem. Namely, 



ffM |t|2. 
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the solution p{x, t) satisfies the wave equation in the whole space, not just 
inside the surface albeit we observe the solution on S only. Since the 
initial function is compactly supported, and since appropriate non-trapping 
condition will be imposed, the energy will be locally decaying (e.g., |26l [38l 
[39]). Thus, the energy of the solution p inside S will be decaying. The rate of 
this decay is different in odd and even dimensions ([SlEHlEn]), exponential in 
the odd and power in the even case. In both cases, the solution is summable 
with respect to time as a function in appropriate functional spaces in the 
domain B bounded by S. We will see that this influences the solvability of 
the problem significantly. 

3.1 A test case: constant sound speed in an odd di- 
mension 

We start with the simplest case when the dimension n is odd and the sound 
speed is constant and equal to 1. We also assume here that the initial function 
f{x) is smooth and compactly supported. Then both the solution t) 
and the boundary data g{x,t) are smooth and, due to Huygens' principle, 
compactly supported in time, when x & B. Here we denote by B the interior 
domain of 5*. Let now E be the operator of harmonic extension of functions 
from S to 5, which, due to standard regularity theorems, is continuous from 
H'^{S) to H^^2(^B) for any s > 0. Then we can rewrite the problem ([2]) in 
terms of the function u = p — Eg: 

' utt - A^u = -E{gtt),t >0,xeB 

^ u{x,0) = f{x)-{Eg\t=o), 
Ut{x, 0) = 

My^t) = o,y e s xR+ 

In particular, if the initial function is supported strictly inside B (which we 
do not need to assume), then g\t=o = 0, and thus u{x, 0) = f{x). 

Let us now denote hy the Dirichlet Laplacian in B, defined in the 
standard way as an unbounded self-adjoint operator in L^{B). 

Theorem 1. The following representation ofu{x,t) holds: 

oo 

u{x,t) = j {~AD)~^^^in{{t-T){-AD)'^)E{gtt){x,T)dT. (5) 
t 



(4) 
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In particular, 

fix) = {Eg\t=o) - I {-ADr^'sm{T{-ADr')E{gu){x,T)dT. (6) 



oo 





Proof. First of all, the integral in ([5]) makes sense, since 

{-ADyhm{{t-T) (-Az))^) 
is an uniformly bounded strongly continuous operator function in LP'{B) and 

a{t) ■.= E{gu)i-,t) 

is a smooth compactly supported function of t with values in L'^{B). Let 
{A^}, be the spectrum and an orthonormal basis of eigenfunctions 

of —Ad, where we assume > 0. We introduce notations 

fk ■= if,i^k),Ukit) := and 0^(0 := 

for the Fourier coefficients of /, u, and a, with respect to the system {i/jk}- 
Here (, ) denotes the scalar product in LP'i^B). 

It is sufficient to prove (l5l) coupled with any eigenfunction ipk'- 

oo 

Uk{t) = j {{-ADy-^sm{{t-T){-AD)^a{r),iJk{x))dT. (7) 
t 

_i 1 
Using self-adjointness of {—Ao) ^ sin ((t — r) (— A/5) ^ ) and the eigenfunc- 
tion property of ipki this simplifies to 



Uk{t) = y sin((t - r)Afe)afe(r)rfr. (8) 

t 

It is easy to check that the right hand side in ([8]) is the unique solution 
vanishing at infinity of the equation 

utt = -Xlu - ak{t). (9) 

However, iQ coincides with the projection of (j4]) onto the eigenfunction i/jk. 
Since the system of eigenfunctions is complete in L^(i?), this proves the 
theorem. □ 
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Let us now try to remove the appearance of the extension operator E. In 
order to do so, we introduce the Fourier coefficients of the data: 



9k{t) =< g, ^ >:= J g{x, t)^^{x)dx, 

s 

where u is the external normal to S. 
Lemma 2. The following equality holds: 

{EgUt) = Xfgkit). (10) 

In particular, 

Mt) = K^g'm- (11) 

Proof. Indeed, {Eg)k = {Eg,ipk) = ^k'^i^g, A^^). Using now Green's for- 
mula and harmonicity of Eg, one gets {Eg)k = A^^ < g, >= X^'^gk(t). □ 

The formulas dHl) and ([10])- ([U]) lead to 



oo 

fk = K^gui^) - J sin {Xkt)g'l{t)dt. (12) 



Now integration by parts in (fT2l) proves the following 

Theorem 3. Function f{x) in ^ can be reconstructed inside B from the 
data g in as the following L'^ (B) -convergent series: 

f{x) = J2fkM^)^ (13) 

k 

where the Fourier coefficients fk can be recovered using one of the following 
formulas: 

oo 

-2^ fn\ \-3 



fk = K 9kiO) - f sin {Xkt)g'l{t)dt, 



oo 

fk = XfgkiO) + Xf J cos {Xkt)g',,{t)dt, or 



fk = ^ I sin {Xkt)gk{t)dt = -A^ ^ / / sin {Xkt)g{x, t)^{x)dxdt, 



dv 
s 



(14) 
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where 




'x 



s 



and u denotes the external normal to S. 

Corollary 4. The statement of the Theorem with the third of the coefficient 
formulas in holds under much milder conditions on the function f{x). 
For instance, it is sufficient to assume that f G //^^^^(M") for some s > |. 

Indeed, one can approximate such a function / in the space //^^^^(R"') 
by functions from C^(M"), thus getting convergence of the data g{x,t) in 
^compi'S X M). Then the third formula in (fT^ extends by continuity. 

3.2 An abstract theorem 

It is not hard to establish an abstract analog of Theorems [H and [3l 

Theorem 5. Let H be a separable Hilbert space and A be a positive self- 
adjoint operator in H with discrete spectrum {A^} and orthonormal basis of 
eigenvectors ipk- Suppose that function u : R"*" H solves the equation 



Utt = —An — a{t) 
and that the following conditions are satisfied: 

1. The function u and its derivative tend to zero when t ^ oo: 



(15) 



lim ||ti(t)||// = lim \\ut{t)\\H = 0. 

t— >oo t— >oo 



2. The function a{t) belongs to L^(R'^,H). 



Then 



1. The following representation of u(t) holds: 



oo 




(16) 



t 



In particular, 



oo 




(17) 
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2. The coefficients Uk of the Fourier series expansion 

k 

can be found as follows: 

oo 

Uk = ~'^k^ J sin(T"Afc)afe(r)dr, 



where ak(t) = {a{t),^pk)H ■ 
Proof. The proof repeats the proofs of Theorems [T] and [31 □ 

3.3 An application to TAT in inhomogeneous media 

Let us consider now the TAT problem ([2]) in the presence of an inhomoge- 
neous background, assuming the conditions imposed on the sound speed c{x) 
in the beginning of Section [31 Using the same harmonic extension trick as in 
the Section [3m we get 

'uu - c^{x)A^u = -E{gu),t >0,xeB 

^ u{x,0) = f{x)-{Eg\t=o), 
Ut{x, 0) = 

^u{y,t) = 0,y e S xR+ 

Consider the Hilbert space H = L'^{B,c~'^{x)dx), i.e., the weighted space 
with the weight c~^(x). On this space, the naturally defined operator 

A = -c\x)A 

in B with zero Dirichlet conditions on S is self-adjoint, positive, and has 
discrete spectrum {A^}(Afc > 0) with eigenfunctions ipk{x) G H. Now, since 
we are in fact dealing with the unobstacled whole space wave propagation 
(the surface S is not truly a boundary, but just an observation surface) and 
since we assumed the sound speed constant at infinity and non-trapping, the 
local energy decay type estimates of [381 ES] (see also [H Theorem 2.104]) 
apply. They, in combination with the standard Sobolev trace theorems, imply 
in particular that the conditions of Theorem [51 both on u{t) := u{-,t) G H 
and a(t) = E{gtt) are satisfied. We thus get the following 
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(18) 



Theorem 6. 1. The function f{x) in /figj). and thus in ^ can he recon- 
structed inside B as follows: 



fix) = {Eg\t=o) - j sm{TA-^)E{gu){x,r)dT. (19) 





2. Function f{x) can be reconstructed inside B from the data g in as 
the following L'^ (B) -convergent series: 



fix) = ^fki^k{x), 



(20) 



where the Fourier coefficients fk can he recovered using one of the fol- 
lowing formulas: 



fk = XfgkiO) - J sin {Xkt)g'at)dt, 



oo 

fk = XfgkiO) + Xf J cos {Xkt)g'k(t)dt, or 



fk = -Xf^^ J sin {Xkt)gk{t)dt = -X,:' J J sm{Xkt)g{x,t)^ix)dxdt, 



s 



(21) 



where 



9k{t) = I g{x,t)^^{x)dx 



and 1/ denotes the external normal to S. 

Proof. The only thing that might have been different is deriving the Fourier 
coefficient representations (j2Tl) by couphng the expression E{gu) with i/jk in 
(fT9l) in the weighted space. This, however, causes no problem, since 



/ E{gtt){x,t)%l)k{x)c ^{x)dx 



XfJ E{gu){x,ty{x)AMx)c-\x)dx 



A^^/ E{gtt){x,t)A'i/;k{x)dx, 

B 



from where the proof proceeds as before. 



□ 
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Uniqueness of reconstruction of f{x) inside B obviously follows from this 
Theorem. It, however, requires some additional arguments, if one wants to 
derive uniqueness of recovery of /(x) outside S. We thus sketch an indepen- 
dent proof of uniqueness of reconstruction of the whole function f{x), which 
does not rely upon the previous theorem. 

Theorem 7. Under the non-trapping conditions formulated above, compactly 
supported function f{x) is uniquelly determined by the data g. (No assump- 
tion of f being supported inside S is imposed.) 

Proof. Suppose that the data g vanishes. Then, inside B, due to local energy 
decay [SIIMIES], u{t), as an //-valued function, is integrable over M. Taking 
Fourier transform in time, we get a continuous //-valued function u{X) that 
satisfies the equation (— -|- A)u{X) = almost everywhere. However, the 
spectrum of A is discrete, and thus u{X) must vanish for all but discrete set 
of values of A. Due to continuity of u{X), this implies that u{X), and thus 
u{t) vanish. This gives vanishing of / inside B. (This could also be derived 
immediately from the preceding theorem). Now, the function u extends to 
a solution of the wave equation outside B with matching Cauchy data on 
S. After Fourier transform in time, we get for each value of A a solution 
u{x,X) of the equation c'^{x)Au + X'^ii = in R", vanishing inside S. Then 
standard uniqueness of continuation theorems (remember that c is smooth 
and non-vanishing) prove that u, and hence u, is zero. □ 

4 Further discussion and remarks 

• The results of this paper are in some sense the realizations of the re- 
marks made in [11] concerning inversion by solving that wave equation 
in reverse time direction and in [2] concerning inversions by eigenfunc- 
tion expansions (improved and implemented in [25]). 

• It is instructive to check how the last of the formulas fll4p works in 
ID in the case of constant speed. Let us assume that B = [0, 1] and 
f{x) is supported in B. Then A^ = kn and ipk{x) = -\/2 sin fcvra;. The 
d'Alambert formula claims that u{x, t) = l{f{x -t) + f{x + t)) and 
thus g{0,t) = \f{t), g{l,t) = |/(1 - t) for t e [0,1] and ^ = for 
t > 1. The normal derivatives of the eigenfunctions are 



(0) 



V2k7v, 



(1) = V2{-lfk7V. 



dv 



dv 
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Now substitution of these expression into the last of the formulas fll4|] 
after a short computation amounts to the standard formula for the 
Fourier coefficients of / with respect to the system ipk- 

• Another instructive computation is to compare with the expansion for- 
mulas from [22] that were obtained under the assumption of a constant 
sound speed. Consider the 3D case. One needs to notice that the data 
used in |25], although denoted by g{t) as ours, is different. In [25], it 
was, up to a factor 47r, the integral over the sphere of radius t, while 
ours corresponds to the solution of the wave equation ([2]). Taking this 



Substitution of this g instead of g into the last formula of ffl^ and one 
integration by parts immediately lead to the formulas (10), (12), and 
(13) in [23]. Note however that, unlike in [25], no knowledge of the 
whole space Green's function is needed. In fact, our formulas hold in 
the case of variable non-trapping sound speed, when one cannot write 
this Green's function explicitly and where the method of [25] fails. 

• Notice that, like in |25|, and unhke [IHl HH 113, the inversion pro- 
cedures derived in this paper do not require the function f{x) to be 
supported inside the observation surface S. Even if the (compact) sup- 
port of / reaches outside S, reconstruction inside S is correct. It is 
known, however, that in this case the backprojection formulas from 
[TUl [TT| [2H W2\ produce incorrect reconstructions inside S. The reason 
is that the formulas from [TOj [TTl [231 112] have the information built 
in about the support of / being inside. In particular, they use, in our 
notations, only the values of t up to the diameter of the domain B. 
It is clear that in this case reconstruction of functions supported not 
necessarily inside S is not unique, since, due to the finite speed of prop- 
agation, information from all support of / might not reach 5* during 
this time. On the other hand, formulas of this paper, as well as of [25] 
use the information for all values of time t. 

• It has been assumed throughout the paper that the initial velocity 
Pt{x,0) in ([2]) is equal to zero. It is clear, however, that this is of no 
importance for our considerations, and one can reconstruct the (possi- 
bly, non-zero) initial velocity as well, using formulas and flTB]) . These 
formulas were derived without using the assumption that pt{x, 0) = 0. 



into account, our data (call it temporarily g), is equal to 
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• We went along with accepted models [3 [TTl EZ] of non-homegenous 
media in TAT that lead to ([2]). In these models, variations of density 
are neglected. However, the abstract Hilbert space formalism that we 
provided, allows one to include also variations of density without much 
of a change in formulas and methods. This just leads to a more general 
operator A in L'^{B). 

• The integral formulas ([5]), f|T6l) . and f|T9|l use functions of the Laplace 
(or a more general) operator in L?{B) with Dirichlet boundary condi- 
tions on 5*. It is clear how to use such formulas in combinations with 
eigenfunction expansions of the operator, as we did in this text. It is, 
however, less clear how such formulas could lead to analytic backpro- 
jection type inversions, such as the ones in [lOl [TTl 1211 US] . A Fourier 
transform consideration that leads to Helmholtz type equations sug- 
gests that there might be a link to [;24], and thus one might hope to 
produce inversion formulas of [25 obtained for the case of S being a 
sphere, as a consequence of our formulas. The following simple com- 
ments might be useful for the further progress in this direction. 

First of all, the gist of the consideration of Section 13.11 is the following 
simple relation (derived readily using Green's formula): 

{ADu)k = {Au)k + J u^dA{x). 
s 



Secondly, reconstruction formulas (!20|) - (pT!) can be combined into the 
following integral formula: 



f{x) = ~J J d,^K{x,y,t)g{y,t)dA{y)dt, (22) 
s 

where the kernel K is 

K{x, y,t) = Y, sm{Xkt)Mx)My)- (23) 

k 

It is readily checked that K satisfies the following initial-boundary value 
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problem: 



'Ktt - c^{y)AyU = 
Kix,y,0) = 0, 
Kt{x,y,0) = S{y 

^K{x,y,t) = 0,y 



eSxR+ 



0,t>0,y e B 



(24) 



In this context, the inversion formula (21) can be formally derived from 
the Green's formula: 



= /bx[o,t][(^« - c\y)AyK)u - K{uu - c\y)Au)]c-\y)dydt 
= Isx[o,T]i^^-y9 - d,^Kg]dA{y)dt + /bx({o}u{t}) t^*^ " Kut]dy, 



by substituting the initial-boundary conditions, letting T ^ oo and 
using the vanishing of the solutions when T — >■ cx). 

It might be possible in the case of S being a sphere to convert the 
inversion formula for / into a backprojection type one. 

• For the constant sound speed case, it has been known for quite some 
time [U [3l [20] that any closed surface S provides uniqueness of re- 
construction. The method of [3] is not applicable, since it relies upon 
constant sound speed. However, spectral methods of [H [20] are appli- 
cable, as it is shown in the proof of Theorem [71 

Theorem [7] guarantees unique recovery of any compactly supported 
function f{x), even if its support is not confined to B. However, mi- 
crolocal arguments show that reconstruction should become unstable 
for some parts of / outside the observation surface S. This is related to 
the existence in this case of the so called "visible" or "audible" zones 
of the wave front set of f{x), as well as those that are not "visible" 
("audible") [271 [Ml [3S]. See [1 [22] for a brief discussion of this 



• Albeit this might not be clear from the text of this paper, what has 
led the authors to this work, was an approach through the range con- 
ditions of the spherical mean operator described in [2J. Indeed, the line 
of thought was that the whole possibility of inversion is based upon 
a very special type of the boundary data g involved. As we have al- 
ready explained, using arbitrary functions as the data g would lead to 
impossibility of reconstruction. Thus, the basis of our approach was 



issue. 
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to use our knowledge of the special features of the data g to derive 
inversions. It can be shown (albeit we do not do this in this text) that 
decay with time of solutions inside B is directly linked (in the case of 
constant sound speed) to the range descriptions of Continuing this 
consideration, one can also obtain some analogs of range descriptions 
for the case of variable sound speed. 
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